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Q ■ Abstract 

CO I A new metaheuristic optimisation algorithm, called Cuckoo Search (CS), 

^^ . was developed recently by Yang and Deb (2009). This paper presents a more 

extensive comparison study using some standard test functions and newly de- 
C) . signed stochastic test functions. We then apply the CS algorithm to solve 

^^ I engineering design optimisation problems, including the design of springs and 

welded beam structures. The optimal solutions obtained by CS are far better 
than the best solutions obtained by an efficient particle swarm optimiser. We 
C^ . will discuss the unique search features used in CS and the implications for fur- 

ther research. 
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Most design optimisation problems in engineering are often higlily nonlinear, involv- 
ing many different design variables under complex constraints. These constraints 
^ ' can be written either as simple bounds such as the ranges of material properties, 

or as nonlinear relationships including maximum stress, maximum deflection, mini- 
mum load capacity, and geometrical configuration. Such nonlinearity often results in 
multimodal response landscape. Subsequently, local search algorithms such as hill- 
climbing and Nelder-Mead downhill simplex methods are not suitable, only global 
algorithms should be used so as to obtain optimal solutions (Deb 1995, Arora 1989, 
Yang 2005, Yang 2008). 

Modern metaheuristic algorithms have been developed with an aim to carry out 
global search, typical examples are genetic algorithms (Glodberg 1989), particle 
swarm optimisation (PSO) (Kennedy and Eberhart 1995, Kennedy et al 2001). The 
efficiency of metaheuristic algorithms can be attributed to the fact that they imitate 
the best features in nature, especially the selection of the fittest in biological systems 
which have evolved by natural selection over millions of years. Two important 
characteristics of metaheuristics are: intensification and diversification (Blum and 
Roll 2003, Gazi and Passino 2004, Yang 2009). Intensification intends to search 



around the current best solutions and select the best candidates or solutions, while 
diversification makes sure that the algorithm can explore the search space more 
efficiently, often by randomization. 

Recently, a new metaheuristic search algorithm, called Cuckoo Search (CS), has 
been developed by Yang and Deb (2009). Preliminary studies show that it is very 
promising and could outperform existing algorithms such as PSO. In this paper, 
we will further study CS and validate it against test functions including stochastic 
test functions. Then, we will apply it to solve design optimisation problems in 
engineering. Finally, we will discuss the unique features of Cuckoo Search and 
propose topics for further studies. 

2 Cuckoo Search 

In order to describe the Cuckoo Search more clearly, let us briefly review the inter- 
esting breed behaviour of certain cuckoo species. Then, we will outline the basic 
ideas and steps of the proposed algorithm. 

2.1 Cuckoo Breeding Behaviour 

Cuckoo are fascinating birds, not only because of the beautiful sounds they can 
make, but also because of their aggressive reproduction strategy. Some species such 
as the ani and Guira cuckoos lay their eggs in communal nests, though they may 
remove others' eggs to increase the hatching probability of their own eggs (Payne et 
al 2005). Quite a number of species engage the obligate brood parasitism by laying 
their eggs in the nests of other host birds (often other species). There are three basic 
types of brood parasitism: intraspecific brood parasitism, cooperative breeding, and 
nest takeover. Some host birds can engage direct conflict with the intruding cuckoos. 
If a host bird discovers the eggs are not its owns, it will either throw these alien eggs 
away or simply abandons its nest and builds a new nest elsewhere. Some cuckoo 
species such as the New World brood-parasitic Tapera have evolved in such a way 
that female parasitic cuckoos are often very specialized in the mimicry in colour and 
pattern of the eggs of a few chosen host species (Payne et al 2005). This reduces the 
probability of their eggs being abandoned and thus increases their reproductivity. 

Furthermore, the timing of egg-laying of some species is also amazing. Parasitic 
cuckoos often choose a nest where the host bird just laid its own eggs. In general, the 
cuckoo eggs hatch slightly earlier than their host eggs. Once the first cuckoo chick 
is hatched, the first instinct action it will take is to evict the host eggs by blindly 
propelling the eggs out of the nest, which increases the cuckoo chick's share of food 
provided by its host bird (Payne et al 2005). Studies also show that a cuckoo chick 
can also mimic the call of host chicks to gain access to more feeding opportunity. 

2.2 Levy Flights 

In nature, animals search for food in a random or quasi-random manner. In general, 
the foraging path of an animal is effectively a random walk because the next move 
is based on the current location/state and the transition probability to the next 
location. Which direction it chooses depends implicitly on a probability which can 
be modelled mathematically. For example, various studies have shown that the flight 



behaviour of many animals and insects has demonstrated the typical characteristics 
of Levy flights (Brown et al 2007, Reynods and Frye 2007, Pavlyukevich 2007). 

A recent study by Reynolds and Frye (2007) shows that fruit flies or Drosophila 
melanogaster, explore their landscape using a series of straight flight paths punc- 
tuated by a sudden 90° turn, leading to a Levy-flight-style intermittent scale-free 
search pattern. Studies on human behaviour such as the Ju/'hoansi hunter-gatherer 
foraging patterns also show the typical feature of Levy flights. Even light can be 
related to Levy flights (Barthelemy et al 2008). Subsequently, such behaviour has 
been applied to optimization and optimal search, and preliminary results show its 
promising capability (Shlesinger 2006, Pavlyukevich 2007). 

2.3 Cuckoo Search 

For simplicity in describing our new Cuckoo Search (Yang and Deb 2009), we now 
use the following three idealized rules: 

• Each cuckoo lays one egg at a time, and dumps it in a randomly chosen nest; 



• 



The best nests with high quality of eggs (solutions) will carry over to the next 
generations; 

• The number of available host nests is fixed, and a host can discover an alien 
egg with a probability pa G [0, 1]. In this case, the host bird can either throw 
the egg away or abandon the nest so as to build a completely new nest in a 
new location. 

For simplicity, this last assumption can be approximated by a fraction pa of the 
n nests being replaced by new nests (with new random solutions at new locations). 
For a maximization problem, the quality or fitness of a solution can simply be 
proportional to the objective function. Other forms of fitness can be defined in a 
similar way to the fitness function in genetic algorithms. 

Based on these three rules, the basic steps of the Cuckoo Search (CS) can be 
summarised as the pseudo code shown in Fig. [TJ 

When generating new solutions a;'*^^' for, say cuckoo i, a Levy flight is performed 

xf^^^ =xf^ + aeLevy{X), (1) 

where a > is the step size which should be related to the scales of the problem 
of interest. In most cases, we can use a = 0(1). The product © means entry-wise 
multiplications. Levy flights essentially provide a random walk while their random 
steps are drawn from a Levy distribution for large steps 

Levy~u = t"^, (1<A<3), (2) 

which has an infinite variance with an infinite mean. Here the consecutive jumps/steps 
of a cuckoo essentially form a random walk process which obeys a power-law step- 
length distribution with a heavy tail. 

It is worth pointing out that, in the real world, if a cuckoo's egg is very similar 
to a host's eggs, then this cuckoo's egg is less likely to be discovered, thus the fitness 
should be related to the difference in solutions. Therefore, it is a good idea to do 
a random walk in a biased way with some random step sizes. A demo version is 
attached in the Appendix (this demo is not published in the actual paper, but as a 
supplement to help readers to implement the cuckoo search correctly). 



Objective function f{x), x = {xi, ...,Xd)'^; 

Initial a population of n host nests Xi (i = 1,2, ...,n); 

while (t <MaxGeneration) or (stop criterion); 

Get a cuckoo (say i) randomly by Levy flights; 

Evaluate its quality /fitness Ff, 

Choose a nest among n (say j) randomly; 

if (F, > F,), 

Replace j by the new solution; 

end 

Abandon a fraction (pa) of worse nests 

[and build new ones at new locations via Levy flights]; 

Keep the best solutions (or nests with quality solutions); 

Rank the solutions and find the current best; 
end while 
Postprocess results and visualisation; 

Figure 1: Cuckoo Search (CS). 

3 Implementation and Validation 

3.1 Validation and Parameter Studies 

It is relatively easy to implement the algorithm, and then we have to benchmark it 
using test functions with analytical or known solutions. There are many benchmark 
test functions and there is no standard list or collection, though extensive descrip- 
tions of various functions do exist in literature (Floudas et al 1999, Hedar 2005, 
Molga and Smutnicki 2005). For example, Michalewicz's test function has many 
local optima 
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fix) = -> sin(xi) sin(^) , (m = 10), (3) 



in the domain < Xj < vr for i = l,2,...,d where d is the number of dimensions. 
The global mimimum /* « -1.801 occurs at (2.20319,1.57049) for d = 2, while 
/^ PS —4.6877 for d = 5. In the 2D case, its 3D landscape is shown Fig. [2l 

The global optimum in 2D can easily be found using Cuckoo Search, and the 
results are shown in Fig. [3] where the final locations of the nests are marked with 
o. Here we have used n = 20 nests, a = 1 and pa = 0.25. From the figure, we can 
see that, as the optimum is approaching, most nests aggregate towards the global 
optimum. In various simulations, we also notice that nests are also distributed at 
different (local) optima in the case of multimodal functions. This means that CS 
can find all the optima simultaneously if the number of nests are much higher than 
the number of local optima. This advantage may become more significant when 
dealing with multimodal and multiobjective optimization problems. 

We have also tried to vary the number of host nests (or the population size n) 
and the probability pa- We have used n = 5, 10, 15, 20, 50, 100, 150, 250, 500 and 
Pa = 0,0.01,0.05,0.1,0.15,0.2,0.25, 0.4,0.5. From our simulations, we found that 
n = 15 to 25 and pa = 0.15 to 0.30 are sufficient for most optimization problems. 
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Figure 2: The landscape of Michaelwicz's 2D function. 

Results and analysis also imply that the convergence rate, to some extent, is not 
sensitive to the parameters used. This means that the fine adjustment of algorithm- 
dependent parameters is not needed for any given problems. Therefore, we will use 
n = 20 and pa = 0.25 in the rest of the simulations, especially for the comparison 
studies presented later. 

3.2 Standard Test Functions 

Various test functions in literature are designed to test the performance of optimiza- 
tion algorithms (Chattopadhyay 1971, Schoen 1993, Shang and Qiu 2006). Any new 
optimization algorithm should also be validated and tested against these benchmark 
functions. In our simulations, we have used the following test functions. 
De Jong's first function is essentially a sphere function 



f{x) = J2xl x,G [-5.12,5.12], 



(4) 



i=l 



whose global minimum f{x^) = occurs at a;* = (0,0, ...,0). Here d is the dimen- 



sion. 



The generalized Rosenbrock's function is given by 

d-l 



fix) = E [(1 - ^^)' + 100(2^^+1 - ^If 



(5) 



i=l 



which has a unique global minimum /^, = at a;* = (1, 1, ..., 1) 




Figure 3: Initial locations of 20 nests in Cuckoo Search, and their final locations are 
marked with o. 



Schwefel's test function is multimodal 

d 



/(^) = H [-a^jsind 



-500 < Xi < 500, 



(6) 



j=i 



whose global minimum /* = -418.9829(i is at x* = 420.9687(z = 1,2, ...,d). 
Ackley's function is also multimodal 



f{x) = -20exp 



0.2, 



1 



Vlk' 



exp[-^cos(27rxi)] + (20 + e), (7) 



'd^ 



i=l 



with the global minimum /* = at a;* = (0, 0, ..., 0) in the range of —32.768 < Xj < 
32.768 where i = 1,2, ...,d. 
Rastrigin's test function 



f{x) = IQd + Y^lxj - 10cos(27rxi)], 



(8) 



j=i 



has a unique global minimum /* = at (0, 0, ..., 0) in a hypercube —5.12 < Xj < 5.12 
where z = 1, 2, ..., d. 

Easom's test function has a sharp tip 

f{x, y) = - cos(x) cos(y) exp[-(x - iif - {y - tt)^], (9) 

in the domain (x, y) G [—100, 100] x [—100, 100]. It has a global minimum of /* = —1 
at (tt, vr) in a very small region. 

Griewangk's test function has many local minima 



/(-) = 77LE-?-n-(^)+i> 



4000^ 



i=l j=l 



^^' 



(10) 



but a unique global mimimum /* = at (0, 0, ..., 0) for all —600 < Xj < 600 where 
i = 1,2, ...,d. 



3.3 Stochastic Test Functions 

Almost all the test functions in literature are deterministic. It is usually more 
difficult for algorithms to deal with stochastic functions. We have designed some 
stochastic test functions for such a purpose. 

The first test function designed by Yang (2010) looks like a standing-wave func- 
tion with a region of defects 

Ed , A2m ^ d 

J you J -y^ -1^"^'/^ -2e-S-i^'(^''-^) ] -Hcos^^i, m = 5, (11) 

4 = 1 

which has many local minima and the unique global minimum /,, = — 1 at a;,,, = 
(vr, TT, ...,7r) for =15 within the domain —20 < Xj < 20 for i = 1,2, ...,d. Here the 
random variables e^ {i = 1, 2, ..., d) are uniformly distributed in (0, 1). For example, 
if all ej are relatively small (say order of 0.05), a snapshot of the landscape in 2D is 
shown in Fig. [H while for higher values such as 0.5 the landscape is different, also 
shown in Fig. HI 

Yang's second test function is also multimodal but it has a singularity 
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\Xi\) 


exp - 


- ) !sin(2;2) 
1=1 



/(^) =\1^ ^i\xi\) exp [ - 2^ sin(2;i )J , (12) 

i=l 1=1 

which has a unique global minimum /* = at a;^, = (0, 0, ...,0) in the domain 
— 27r < Xi < 2tt where i = l,2,...,d (Yang 2010). This function is singular at the 
optimum (0, ..., 0). Similarly, e^ should be drawn from a uniform distribution in [0, 1] 
or Unif[0,l]. In fact, using the same methodology, we can turn many determistic 
functions into stochastic test functions. For example, we can extend Robsenbrock's 
function as the following stochastic function 

d-l 

fix) = E [(1 - ^*)' + 100e,(xi+i - xf)2] , (13) 

i=l 

where e^ should be drawn from Unif[0,l]. Similarly, we can also extend De Jong's 
function into its corresponding stochastic form 

d 
fix) = Y,eixl (14) 

i=l 

which still has the same global minimum /* = at (0, 0, ..., 0), despite its stochastic 
nature due to the factor e^. For stochastic functions, most deterministic algorithms 
such as hill climbing and Nelder-Mead downhill simplex method would simply fail. 
However, we can see later that most metaheuristic algorithms such as PSO and CS 
are still robust. 

3.4 Simulations and Comparison 

Recent studies indicate that PSO can outperform genetic algorithms (GA) and other 
conventional algorithms (Goldberg 1989, Kennedy et al 2001, Yang 2008). This can 
be attributed partly to the broadcasting ability of the current best estimates, po- 
tentially leading to a better and quicker convergence rate towards the optimality. A 
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Figure 4: Landscape of stochastic function (jlip for small e (left) and large e (right). 
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Figure 5: The initial location of 20 nests (left) for function (|lip and their final 
locations after 15 iterations (right). 



general framework for evaluating statistical performance of evolutionary algorithms 
has been discussed in detail by Shilane et al (2008). 

Now we can compare the Cuckoo Search with PSO and genetic algorithms for 
various test functions. After implementing these algorithms using Matlab, we have 
carried out extensive simulations and each algorithm has been run at least 100 
times so as to carry out meaningful statistical analysis. The algorithms stop when 
the variations of function values are less than a given tolerance e < 10~^. The results 
are summarised in Table 1 where the numbers are in the format: average number of 
evaluations it one standard deviation (success rate), so 3321ib519(100%) means that 
the average number (mean) of function evaluations is 3321 with a standard deviation 
of 519. The success rate of finding the global optima for this algorithm is 100%. 
The functions used in the Table are (1) Michaelwicz (d = 16), (2) Rosenrbrock 
{d = 16), (3) De Jong {d = 32), (4) Schwefel {d = 32), (5) Ackley {d = 128), 
(6) Rastrigin, (7) Easom, (8) Griewangk, (9) Yang's first stochastic function, (10) 
Yang's second stochastic function, (11) Generalised Robsenbrock's function with 
stochastic components, and (12) De Jong's stochastic function. 

We can see that CS is much more efficient in finding the global optima with 



Table 1: Comparison of CS with genetic algorithms and particle swarm optimisation 



Functions 



GA 



PSO 



CS 



(1) 89325 ± 7914(95%) 

(2) 55723 ±8901(90%) 

(3) 15232 ± 1270(100%) 

(4) 23790 ± 6523(95%) 

(5) 32720 ± 3327(90%) 

(6) 110523 ± 5199(77%) 

(7) 19239 ± 3307(92%) 

(8) 70925 ± 7652(90%) 

(9) 79025 ± 6312(49%) 

(10) 35072 ± 3730(54%) 

(11) 63268 ±5091(40%) 

(12) 24164 ± 4923(68%) 



6922 ± 537(98%) 
32756 ± 5325(98%) 
10079 ± 970(100%) 
92411 ± 1163(97%) 
23407 ± 4325(92%) 
79491 ± 3715(90%) 
17273 ± 2929(90%) 
55970 ± 4223(92%) 
34056 ± 4470(90%) 
22360 ± 2649(92%) 
49152 ± 6505(89%) 
11780 ± 4912(94%) 



3221 ± 519(100%) 
5923 ± 1937(100%) 
3015 ± 540(100%) 
4710 ± 592(100%) 
4936 ± 903(100%) 
10354 ± 3755(100%) 
6751 ± 1902(100%) 
10912 ±4050(100%) 
11254 ± 2733(99%) 
8669 ± 3480(98%) 
10564 ± 4297(99%) 
7723 ± 2504(100%) 



higher success rates. Each function evaluation is virtually instantaneous on a modern 
personal computer. For example, the computing time for 10,000 evaluations on a 
3GHz desktop is about 5 seconds. In addition, for stochastic functions, genetic 
algorithms do not perform well, while PSO is better. However, CS is far more 
promising. 

4 Engineering Design 

Design optimisation is an integrated part of designing any new products in engineer- 
ing and industry. Most design problems are complex and multiobjective, sometimes 
even the optimal solutions of interest do not exist. In order to see how the CS 
algorithm may perform, we now use two standard but well-known test problems. 

4.1 Spring Design Optimisation 

Tensional and/or compressional springs are used widely in engineering. A standard 
spring design problem has three design variables: the wire diameter w, the mean 
coil diameter d, and the length (or number of coils) L. 

The objective is to minimise the weight of the spring, subject to various con- 
straints such as maximum shear stress, minimum deflection, and geometrical limits. 
For detailed description, please refer to earlier studies (Belegundu 1982, Arora 1989, 
Cagnina et al 2008). This problem can be written compactly as 



Minimise f{x) = (L ± 2)w^d, 



(15) 



subject to 



9iix) = 1 
92ix) = 1 



71785«)4 — "' 
140. 45«) ^ n 



IPT" 



53(2;) 



2(w+d) 



(16) 



1 <0, 



d{4d-w) 



+ 



1 <0, 



9A^) " «,3(12566d-«i) "•" 5108«)2 

with the following limits 

0.05 <u;< 2.0, 0.25 <(i< 1.3, 2.0 < L < 15.0. 



(17) 



Using Cuckoo Search, we have obtained the same or slightly better solutions 
than the best solution obtained by Cagnina et al (2008) 

/* = 0.012665 at (0.051690,0.356750,11.287126), (18) 

but cuckoo search uses significantly fewer evaluations. 

4.2 Welded Beam Design 

The so-called welded beam design is another standard test problem for constrained 
design optimisation (Ragsdell and Phillips 1976, Cagnina et al 2008). The problem 
has four design variables: the width w and length L of the welded area, the depth 
h and thickness h of the main beam. The objective is to minimise the overall 
fabrication cost, under the appropriate constraints of shear stress r, bending stress 
a, buckling load P and maximum end deflection 6. 
The problem can be written as 



minimise f{x) = 1.10A71w'^L + 0.04811d/i(14.0 + L), 



(19) 



subject to 



9i{x) 


= w-h<0, 


92{x) 


= Six)- 0.25 < 0, 


93{x) 


= r(a;)-13,600 < 0, 


94ix) 


= c7(a;)-30,000<0. 


9b{x) 


= 0.10471^2 + 0.0481 l/i(i(14 + L 


9&{x) 


= 0.125 -w < 0, 


97{x) 


= 6000 - P{x) < 0, 



(20) 



5.0 < 0, 



10 



where 



aH = 5^, Q = 6000(14 + f), 

r _ 65,856 a — QD 

" ~ 30,000M^ ' P - J ^ (21) 



^ 6000 ^(^-^ _ ^ /^2 _L M^ _L fl2 

V2wL' 



t{x) = ^a^ + ^ + /3^ 



P = 0.61423 X 106 ^(1 - ^^^J^; 



The simple hmits or bounds are 0.1 < L, d < 10 and 0.1 < w,h < 2.0. 
Using our Cuckoo Search, we have the following optimal solution 

x^ = {w,L,d, h) 

= (0.205729639786079, 3.470488665627977, 9.036623910357633, 0.205729639786079), 

(22) 
with 

f{x*)r,nn = 1.724852308597361. (23) 

This solution is exactly the same as the solution obtained by Cagnina et al (2008) 

/* = 1.724852 at (0.205730,3.470489,9.036624,0.205729). (24) 

We have seen that, for both test problems, CS has found the optimal solutions which 
are either better than or the same as the solutions found so far in literature. 

5 Discussions and Conclusions 

From the comparison study of the performance of CS with GAs and PSO, we know 
that our new Cuckoo Search in combination with Levy flights is very efficient and 
proves to be superior for almost all the test problems. This is partly due to the fact 
that there are fewer parameters to be fine-tuned in CS than in PSO and genetic algo- 
rithms. In fact, apart from the population size n, there is essentially one parameter 
Pa- If we look at the CS algorithm carefully, there are essentially three compo- 
nents: selection of the best, exploitation by local random walk, and exploration by 
randomization via Levy flights globally. 

The selection of the best by keeping the best nests or solutions is equivalent 
to some form of elitism commonly used in genetic algorithms, which ensures the 
best solution is passed onto the next iteration and there is no risk that the best 
solutions are cast out of the population. The exploitation around the best solutions 
is performed by using a local random walk 

£C*+^ = x^ + aet. (25) 

If St obeys a Gaussian distribution, this becomes a standard random walk indeed. 
This is equivalent to the crucial step in pitch adjustment in Harmony Search (Geem 
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et al 2001, Yang 2009). If St is drawn from a Levy distribution, the step of move 
is larger, and could be potentially more efficient. However, if the step is too large, 
there is risk that the move is too far away. Fortunately, the elitism by keeping the 
best solutions makes sure that the exploitation moves are within the neighbourhood 
of the best solutions locally. 

On the other hand, in order to sample the search space effectively so that new 
solutions to be generated are diverse enough, the exploration step is carried out in 
terms of Levy flights. In contrast, most metaheuristic algorithms use either uniform 
distributions or Gaussian to generate new explorative moves (Geem et al 2001, Blum 
and Rilo 2003). If the search space is large. Levy flights are usually more efficient. 
A good combination of the above three components can thus lead to an efficient 
algorithm such as Cuckoo Search. 

Furthermore, our simulations also indicate that the convergence rate is insensi- 
tive to the algorithm-dependent parameters such as pa- This also means that we 
do not have to fine tune these parameters for a specific problem. Subsequently, CS 
is more generic and robust for many optimisation problems, comparing with other 
metaheuristic algorithms. 

This potentially powerful optimisation strategy can easily be extended to study 
multiobjecitve optimization applications with various constraints, including NP- 
hard problems. Further studies can focus on the sensitivity and parameter studies 
and their possible relationships with the convergence rate of the algorithm. In 
addition, hybridization with other popular algorithms such as PSO will also be 
potentially fruitful. More importantly, as for most metaheuristic algorithms, math- 
ematical analysis of the algorithm structures is highly needed. At the moment, no 
such framework exists for analyzing metaheuristics in general. Any progress in this 
area will potentially provide new insight into the understanding of how and why 
metaheuristic algorithms work. 
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Appendix: Demo Implementation 

I 

7o Cuckoo Search (CS) algorithm by Xin-She Yang and Suash Deb 7, 
7o Programmed by Xin-She Yang at Cambridge University % 
7. Programming dates: Nov 2008 to June 2009 7, 
7. Last revised: Dec 2009 (simplified version for demo only) 7. 
% 

7« Papers — Citation Details: 

7o 1) X.-S. Yang, S. Deb, Cuckoo search via Levy flights, 

7o in: Proc. of World Congress on Nature & Biologically Inspired 

7. Computing (NaBIC 2009), December 2009, India, 

7 IEEE Publications, USA, pp. 210-214 (2009). 

7 http : //arxiv . org/PS_cache/arxiv/pdf /1003/1003 . 1594vl . pdf 

7 2) X.-S. Yang, S. Deb, Engineering optimization by cuckoo search, 

7 Int. J. Mathematical Modelling and Numerical Optimisation, 

7 Vol. 1, No. 4, 330-343 (2010). 

7 http : //arxiv . org/PS_cache/arxiv/pdf /1005/1005 . 2908v2 . pdf 

I 1 

°L This demo program only implements a standard version of 7. 

7 Cuckoo Search (CS) , as the Levy flights and generation of 7 

7 new solutions may use slightly different methods. °L 

7 The pseudo code was given sequentially (select a cuckoo etc) , 7 

7 but the implementation here uses Matlab's vector capability, 7 

7 which results in neater/better codes and shorter running time. 7 

7 This implementation is different and more efficient than the 7o 
7 the demo code provided in the book by 

7 "Yang X. S., Nature-Inspired Metaheuristic Algoirthms, 'L 

7 2nd Edition, Luniver Press, (2010). " 7. 
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I I 

X =============================================================== I 

7. Notes: '/, 
7o Different implementations may lead to slightly different 7, 
7o behavour and/or results, but there is nothing wrong with it, 7o 
7. as this is the nature of random walks and all metaheuristics . 7o 
I 

function [bestnest ,fmin] =cuckoo_search(n) 

if nargin<l, 

7o Number of nests (or different solutions) 

n=25; 

end 

7o Discovery rate of alien eggs/solutions 
pa=0.25; 

7i,7i, Change this if you want to get better results 

7b Tolerance 

Tol=1.0e-5; 

7o7o Simple bounds of the search domain 

7o Lower bounds 

nd=15; 

Lb=-5*ones(l ,nd) ; 

7o Upper bounds 

Ub=5*ones(l,nd) ; 

7« Random initial solutions 

for i=l:n, 

nest(i, : )=Lb+(Ub-Lb) . *rand(size(Lb) ) ; 

end 

7o Get the current best 

f itness=10"10*ones(n, 1) ; 

[f min , bestnest , nest , fitness] =get_best_nest (nest ,nest , fitness) ; 

N_iter=0; 

7„7„ Starting iterations 

while (fmin>Tol) , 

7o Generate new solutions (but keep the current best) 
new_nest=get_cuckoos (nest , bestnest , Lb, Ub) ; 

[f new , best , nest , fitness] =get_best_nest (nest , new_nest , fitness) ; 
7o Update the counter 

N_iter=N_iter+n; 
7o Discovery and randomization 

new_nest=empty_nests(nest ,Lb,Ub,pa) ; 

7o Evaluate this set of solutions 

[fnew,best ,nest , fitness] =get_best_nest (nest ,new_nest , fitness) ; 
7a Update the counter again 

N_iter=N_iter+n; 
7o Find the best objective so far 
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if fnew<fmin, 
f min=f new ; 
bestnest=best ; 
end 
end VL End of iterations 

VL Post-optimization processing 

°/°L Display all the nests 

dispCstrcat ( 'Total number of iterations=' ,num2str (N_iter)) ) ; 

fmin 

bestnest 

°/°L All subfunctions are list below 

°/°L Get cuckoos by ramdom walk 

function nest=get_cuckoos(nest ,best ,Lb,Ub) 

"/o Levy flights 

n=size(nest , 1) ; 

"/o Levy exponent and coefficient 

"/o For details, see equation (2.21), Page 16 (chapter 2) of the book 

°/o X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Luniver Press, (2010) 

beta=3/2; 

sigma= (gamma (1+beta) *sin(pi*beta/2) /(gamma( (1+beta) /2) *beta*2~ ( (beta-1) /2) ) ) ~ (1/beta) ; 

for j=l:n, 

s=nest (j , : ) ; 

7o This is a simple way of implementing Levy flights 

7o For standard random walks, use step=l; 

7o7o Levy flights by Mantegna's algorithm 

u=randn(size(s) )*sigma; 

v=randn(size(s) ) ; 

step=u. /abs(v) . ~ (1/beta) ; 

7o In the next equation, the difference factor (s-best) means that 

7o when the solution is the best solution, it remains unchanged. 

stepsize=0 . 01*step . * (s-best) ; 

7o Here the factor 0.01 comes from the fact that L/lOO should the typical 

7o step size of walks/flights where L is the typical lenghtscale; 

7o otherwise. Levy flights may become too aggresive/ef f icient , 

7o which makes new solutions (even) jump out side of the design domain 

7o (and thus wasting evaluations) . 

7o Now the actual random walks or flights 

s=s+stepsize . *randn(size(s) ) ; 
°L Apply simple bounds/limits 
nest (j , : )=simplebounds(s ,Lb,Ub) ; 
end 

7o7o Find the current best nest 

function [fmin, best ,nest , fitness] =get_best_nest (nest ,newnest , fitness) 

7o Evaluating all new solutions 

for j=l : size(nest , 1) , 

f new=f ob j (newnest ( j , : ) ) ; 
if f new<=f itness(j) , 
f itness ( j ) =f new ; 
nest ( j , : ) =newnest ( j , : ) ; 
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end 
end 

°L Find the current best 
[fmin,K] =min(f itness) ; 
best=nest (K, : ) ; 

°/o7o Replace some nests by constructing new solutions/nests 

function new_nest=empty_nests(nest ,Lb,Ub,pa) 

"/o A fraction of worse nests are discovered with a probability pa 

n=size(nest , 1) ; 

°/o Discovered or not — a status vector 

K=rand(size(nest) )>pa; 

°h In the real world, if a cuckoo's egg is very similar to a host's eggs, then 
*/ this cuckoo's egg is less likely to be discovered, thus the fitness should 
% be related to the difference in solutions. Therefore, it is a good idea 
"/o to do a random walk in a biased way with some random step sizes. 
°/o°/o New solution by biased/selective random walks 
stepsize=rand* (nest (randperm(n) , : ) -nest (randperm(n) , : ) ) ; 
new_nest=nest+stepsize . *K ; 

°/o Application of simple constraints 
function s=simplebounds(s,Lb,Ub) 

7o Apply the lower bound 

ns_tmp=s ; 

I=ns_tmp<Lb; 

ns_tmp(I)=Lb(I) ; 

7o Apply the upper bounds 

J=ns_tmp>Ub; 

ns_tmp(J)=Ub(J) ; 

7o Update this new move 

s=ns_tmp; 

7o7o You can replace the following by your own functions 

7o A d-dimensional objective function 

function z=f obj (u) 

7o7o d-dimensional sphere function sum_j=l~d (u_j-l)~2. 

7o with a minimum at (1,1, . . . . , 1) ; 

z=sum((u-l) . ~2) ; 
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